#include<stdio.h>


#include "Math/LB_Math.h"
#include "Matrix/LB_Matrix.h"
#include "Algorithm/subdivision_of_polygon.h"
#include "tools_node.h"
#include <Arcroll.h>
#include <directories_and_files.h>
#include<Mesh_IO/Mesh_IO.h>
#include<Mesh/libcell_iterator.h>
#include <time.h>
#include<libcell_macros.h>
#include<Algorithm/threeD_elementary_geometry_tool_library.h>
#include <test_two_method.h>
int global_size=1200;

static double* make_dense_remove_points(double* points, double max_level,double* boundings);

static inline  int point_is_in_triangle(double *p1,double * p2,double* p3,double *p);



double* get_bounding_box_of_one_triangle(double* p1,double*  p2,double* p3)
{
    double* ps[3]={p1,p2,p3};
    double * re=(double*)malloc(sizeof(double)*4);
    re[0]=p1[0];re[1]=p1[1];re[2]=p1[0];re[3]=p1[1];
    for(int i=0;i<3;i++)
    {
        for(int j=0;j<2;j++)
        {
            re[j]= (ps[i][j]<re[j]?(ps[i][j]):re[j]);
            re[j+2]= (ps[i][j]>re[j+2]?(ps[i][j]):re[j+2]);
        }
    }
    return re;
}



static inline  int point_is_in_triangle(double *p1,double * p2,double* p3,double *p)
{
    double t1[3]={p[0]-p1[0],p[1]-p1[1],p[2]-p1[2]};
    double t2[3]={p[0]-p2[0],p[1]-p2[1],p[2]-p2[2]};
    double t3[3]={p[0]-p3[0],p[1]-p3[1],p[2]-p3[2]};
    double* n1=out_product(t1,t2);
    double* n2=out_product(t2,t3);
    double* n3=out_product(t3,t1);
    int re=0;
    if(inner_product(n1,n2,3)>0&&inner_product(n1,n3,3)>0 )
    {
        re= 1;
    }
    SAFE_FREE(n1);SAFE_FREE(n2);SAFE_FREE(n3);
    return re;
}

double* get_bounding_box_of_one_polygon(double**ps,int len)
{
    double  * re=(double*)malloc(sizeof(double)*4);
    re[0]=ps[1][0];re[1]=ps[1][1];re[2]=ps[1][0];re[3]=ps[1][1];
    for(int i=0;i<len;i++)
    {
        for(int j=0;j<2;j++)
        {
            re[j]= (ps[i][j]<re[j]?(ps[i][j]):re[j]);
            re[j+2]= (ps[i][j]>re[j+2]?(ps[i][j]):re[j+2]);
        }
    }
    return re;
}
double**  pre_get_pixel_data(double *boundings,double**ps,int len,int ** subs)
{
    //int size=150;
    int size=global_size;
    double** re=(double**)malloc(sizeof(double*)*size);
    for(int i=0;i<size;i++)
    {
        re[i]=(double*)malloc(sizeof(double)*size);
        memset(re[i],0,sizeof(double)*size);
    }   
    double stepx=(boundings[2]-boundings[0])/((double)size);
    double stepy=(boundings[3]-boundings[1])/((double)size);

    double area=stepx*stepy,sum=0,p[3]={0,0,0}; 
    double* boundings1=NULL;
    int minx=0,miny=0,maxx=0, maxy=0;   
    for(int i=0;i<len-2;i++ )
    {
        boundings1=get_bounding_box_of_one_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]]);
        minx=(int)((boundings1[0]-boundings[0])/(stepx));
        miny=(int)((boundings1[1]-boundings[1])/(stepy));
        maxx=(int)((boundings1[2]-boundings[0])/(stepx))+1;
        maxy=(int)((boundings1[3]-boundings[1])/(stepy))+1;
        minx=((minx<0)?0:(minx));
        miny=((miny<0)?0:(miny));
        maxx=((maxx>size)?(size):maxx);
        maxy=((maxy>size)?(size):maxy);
        
        for(int k=minx;k<maxx;k++)
        {
            p[0]=boundings[0]+stepx*(k+0.5);
            for(int j=miny;j<maxy;j++)
            {
                p[1]=boundings[1]+stepy*(j+0.5);
                if(re[k][j]!=-1&&point_is_in_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]],p))
                {
                    re[k][j]=-1;
                }
            }
        }
        free(boundings1);
    }
    return re;
}

double** get_pixel_datan(double*boundings,double**data,Node* n)
{
    //int size=150;
    int size=global_size;
    double** re=(double**)malloc(sizeof(double*)*size);
    for(int i=0;i<size;i++)
    {
        re[i]=(double*)malloc(sizeof(double)*size);
        memset(re[i],0,sizeof(double)*size);
    }
    double stepx=(boundings[2]-boundings[0])/((double)size);
    double stepy=(boundings[3]-boundings[1])/((double)size);

    double area=stepx*stepy,sum=0,p[3]={0,0,0},np[3]={0,0,0},xs[3]={0,0,0}; 
    double* boundings1=NULL;
    double temp=0;
    int minx=0,miny=0,maxx=0, maxy=0;   

    LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,3,3); 
    double*mat_data=(double*)(mat->data);
    mat_data[0*3+2]=0; mat_data[1*3+2]=0; mat_data[2*3+2]=1.0;
    for(Node*nit=n;nit!=NULL;nit=(Node*)(nit->Next) )
    {
        double* tri_c=(double*)(nit->value);
        for(int j=0;j<3;j++)
        {
            mat_data[j*3+0]=tri_c[1*3+j]-tri_c[0*3+j]; // c->vertices[1]->point[j]-c->vertices[0]->point[j];
            mat_data[j*3+1]= tri_c[2*3+j]-tri_c[0*3+j];    //  c->vertices[2]->point[j]-c->vertices[0]->point[j];
        }   
        LB_Matrix* inv_mat=mat->inverse(mat);
        if(inv_mat==NULL)
        {
            continue;
        }
        double *inv_mat_data=(double*)(inv_mat->data);
        boundings1=get_bounding_box_of_one_triangle(&tri_c[0*3],&tri_c[1*3],&tri_c[2*3]);
        minx=(int)((boundings1[0]-boundings[0])/(stepx));
        miny=(int)((boundings1[1]-boundings[1])/(stepy));
        maxx=(int)((boundings1[2]-boundings[0])/(stepx))+1;
        maxy=(int)((boundings1[3]-boundings[1])/(stepy))+1;
        minx=((minx<0)?0:(minx));
        miny=((miny<0)?0:(miny));
        maxx=((maxx>size)?(size):maxx);
        maxy=((maxy>size)?(size):maxy);
        np[2]=p[2]-tri_c[0*3+2];
        for(int i=minx;i<maxx;i++)
        {
            p[0]=boundings[0]+stepx*(i+0.5);
            np[0]=p[0]-tri_c[0*3+0];
            for(int j=miny;j<maxy;j++)
            {
                if(data[i][j]==0)
                {
                    continue;
                }
                p[1]=boundings[1]+stepy*(j+0.5);
                np[1]=p[1]-tri_c[0*3+1];
                //此处应该是对单个三角形
                xs[0]=xs[1]=xs[2]=0;
                for(int k=0;k<3;k++)
                {
                    xs[0]+=np[k]*inv_mat_data[0*3+k];
                    xs[1]+=np[k]*inv_mat_data[1*3+k];
                }
                if(xs[0]>=0&&xs[1]>=0&&xs[0]+xs[1]<=1)
                {
                    temp=xs[0]*tri_c[1*3+2]+ xs[1]*tri_c[2*3+2]+(1-xs[0]-xs[1])*tri_c[0*3+2];
                    
                    if(data[i][j]==-1)
                    {                       
                        re[i][j]=temp;
                        data[i][j]=-2;
                    }
                    else if(temp>re[i][j])
                    {
                        re[i][j]=temp;
                    }
                }
            }
        }
        lb_matrix_free(inv_mat);
        free(boundings1);
    }
    lb_matrix_free(mat);
    for(int i=0;i<size;i++)
    {
        for(int j=0;j<size;j++)
        {
            data[i][j]==-2?(data[i][j]=-1):(data[i][j]=data[i][j]);
        }
    }
    return re;
}
Node* get_some_triangles_of_intersections_from_mesh(double* boundings,double* m,int len_m)
{
    Node* re=NULL;
    int num=len_m/9;
    for(int i=0;i<num;i++)
    {
        for(int j=0;j<3;j++)
        {
            if( boundings[0]<=m[i*9+j*3+0]&& m[i*9+j*3+0]<=boundings[2]&&  boundings[1]<=m[i*9+j*3+1]&&m[i*9+j*3+1]<=boundings[3] )
            {
                re=node_overlying(re,&m[i*9]);
                break;
            }   
        }   
    }   
    return re;
}


static double*  my_test3( double* m1,int len_m1,double**ps ,int len )
{
    //printf("begin\n");
    int**subs=(int**)malloc(sizeof(int*)*(len-2));
    for(int i=0;i<len-2;i++)
    {
        subs[i]=(int*)malloc(sizeof(int)*3);
        memset(subs[i],0,sizeof(int)*3);
    }
    double max_level=ps[0][2];
    for(int i=0;i<len;i++)
    {
        if(ps[i][2]>max_level)
        {
            max_level=ps[i][2];
        }
    }
    for(int i=0;i<len;i++)
    {
        ps[i][2]=0;
    }
    printf("max_level:%lf\n",max_level);

    optimization_subdivision_of_polygon(ps,len,subs);
    double* boundings =get_bounding_box_of_one_polygon(ps,len);
    Node*n1=get_some_triangles_of_intersections_from_mesh(boundings,m1,len_m1);
    double**pixel_data1= pre_get_pixel_data(boundings,ps,len,subs);
    double**int_pixel_data1= get_pixel_datan(boundings,pixel_data1,n1);
   // int size=150;
    int size=global_size;
    double * re=(double*)malloc(sizeof(double)*(3*size*size+2));
    memset(re,0,sizeof(double)*(3*size*size+2));
    
    double stepx=(boundings[2]-boundings[0])/((double)size);
    double stepy=(boundings[3]-boundings[1])/((double)size);
    double area=stepx*stepy,temp=0,sum=0,sum1=0;int temp_num=0;
    
    for(int i=0;i<size;i++)
    {
        temp=(i+0.5)*stepx+boundings[0];
        for(int j=0;j<size;j++)
        {
            if(pixel_data1[i][j]==-1)
            {
                re[temp_num*3+0 +2]=temp;
                re[temp_num*3+1 +2]=(j+0.5)*stepy+boundings[1];
                re[temp_num*3+2+2]=int_pixel_data1[i][j];
                if(int_pixel_data1[i][j]<max_level)
                {
                    sum+=(int_pixel_data1[i][j]-max_level);
                }
                sum1+=(int_pixel_data1[i][j]-max_level); 
                // printf("%lf   ", (int_pixel_data1[i][j]-max_level));
                // sum+=(int_pixel_data1[i][j]-max_level);
                temp_num++;
            }   
        }
        SAFE_FREE(pixel_data1[i]);
        SAFE_FREE(int_pixel_data1[i]);
    }

    re[0]=area*sum;
    re[1]=temp_num*3;
    printf("结果1:%lf 结果2:%lf \n",area*sum,area*sum1);
 
   // printf(" re :%lf %lf\n",re[0],re[1]);
    double* re1=make_dense_remove_points(re,max_level,boundings);
    free(re);
    SAFE_FREE(int_pixel_data1);
    SAFE_FREE(pixel_data1);
    free_node(n1);
    for(int i=0;i<len-2;i++)
    {
        free(subs[i]);
    }
    free(subs);
   // printf(" ddarea %lf %lf end\n",area*temp_num, compute_area_of_triangle(ps[0],ps[1],ps[2])+compute_area_of_triangle(ps[2],ps[3],ps[0]));
    return re1;
}
static double* make_dense_remove_points(double* points, double max_level,double* boundings)
{
    //int size=150;
    int size=global_size;

    double* re=(double*)malloc(sizeof(double)*(3*size*size*4+2));
    re[0]=points[0];
    int len=points[1]/3;
    int sum=0;
    double stepx=(boundings[2]-boundings[0])/((double)size);
    double stepy=(boundings[3]-boundings[1])/((double)size);
    stepx/=2.0;
    stepy/=2.0; 
    for(int i=0;i<len;i++)
    {
        //printf("fou\n");
        if(points[i*3+2+2]<max_level)
        {
            //printf("shi\n");
            re[sum*3+0+2]=points[i*3+0+2];re[sum*3+1+2]=points[i*3+1+2];re[sum*3+2+2]=points[i*3+2+2];
            sum++;
            re[sum*3+0+2]=points[i*3+0+2]+stepx;re[sum*3+1+2]=points[i*3+1+2];re[sum*3+2+2]=points[i*3+2+2];
            sum++;
            re[sum*3+0+2]=points[i*3+0+2];re[sum*3+1+2]=points[i*3+1+2]+stepy;re[sum*3+2+2]=points[i*3+2+2];
            sum++;
            re[sum*3+0+2]=points[i*3+0+2]+stepx;re[sum*3+1+2]=points[i*3+1+2]+stepy;re[sum*3+2+2]=points[i*3+2+2];
            sum++;
        } 
    } 
    re[1]=sum*3;
    return re;
}

double*   compute_volume_difference_of_two_mesh( double* m1,int len_m1,double*ps ,int len )
{
    clock_t t = clock();
 
    int len1=len/3;
    double** ps1=(double**)malloc(sizeof(double*)*len1);
    for(int i=0;i<len1;i++)
    {
        ps1[i]=(double*)malloc(sizeof(double)*3);
        ps1[i][0]=ps[i*3+0];ps1[i][1]=ps[i*3+1];ps1[i][2]=ps[i*3+2];
    }
    double * re=my_test3(m1,len_m1,ps1,len1);


    for(int i=0;i<len1;i++)
    {
        SAFE_FREE(ps1[i]);
    }
    SAFE_FREE(ps1);
    re[0]=fabs(re[0]);
    //printf("re :%lf %lf\n",re[0],re[1]);
    printf("结果:%lf\n",re[0]);

    printf("运行时间%.4f秒:%d\n", (double)(clock() - t) / CLOCKS_PER_SEC,0); 

    return re;
}

double*  compute_volume_difference_of_two_mesh1( double* m1,int len_m1,double*ps ,int len )
{
    clock_t t = clock();

    int len1=len/3;
    double** ps1=(double**)malloc(sizeof(double*)*len1);
    for(int i=0;i<len1;i++)
    {
        ps1[i]=(double*)malloc(sizeof(double)*3);
        ps1[i][0]=ps[i*3+0];ps1[i][1]=-ps[i*3+2];ps1[i][2]=ps[i*3+1];
    }
    double temp_y=0,temp_z=0;
    for(int i=0;i<len_m1/3;i++)
    {
        temp_y=-m1[i*3+2];
        temp_z=m1[i*3+1];
        m1[i*3+1]=temp_y;
        m1[i*3+2]=temp_z;       

    }
    double * re=my_test3(m1,len_m1,ps1,len1);   
    //my_test3(m1,len_m1,ps1,len1,re);  
    
    for(int i=0;i<re[1]/3;i++)
    {
        temp_y= re[i*3+2+2];
        temp_z=-re[i*3+1+2];
        re[i*3+1+2]=temp_y;
        re[i*3+2+2]=temp_z; 
    }   
    re[0]=fabs(re[0]);
    for(int i=0;i<len1;i++)
    {
        SAFE_FREE(ps1[i]);
    }
    SAFE_FREE(ps1);
    printf("结果:%lf\n",re[0]);
 
    printf("运行时间%.4f秒:%d\n", (double)(clock() - t) / CLOCKS_PER_SEC,0); 

    return re;
}
static void test_show_points_set(Viewer_World*vw,double* points,int len)
{
    Node* n=vw->create_something(vw,"Points");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Points*vp= (Viewer_Points*)(vs->evolution);
    vp->pointsize=10.0;
    vp->Data_rows=len/3;
    vp->Data=(float*)malloc(sizeof(float)*len);
    vp->set_color(vp,0.9,0.9,0.9,1.0);
    for(int i=0;i<len;i++)
    {
        vp->Data[i]=points[i];
    }
    
    free_node(n);
}

static  void test_show_mesh_triangles(Viewer_World* vw,Node* node)
{
    Node* n=vw->create_something(vw,"Faces");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
     
    vf->Data_rows=node_size(node)*3;
    vf->Data_index_rows=node_size(node);
    vf->Data= (float*)malloc(sizeof(float)*vf->Data_rows*3);
    vf->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*vf->Data_index_rows*4);
    int i=0;
    for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        double* tri_c=(double*)(nit->value);
        for(int j=0;j<3;j++)
        {
            vf->Data[i*9+j*3+0]=tri_c[j*3+0];
            vf->Data[i*9+j*3+1]=tri_c[j*3+1];
            vf->Data[i*9+j*3+2]=tri_c[j*3+2];
        }   
        vf->Data_index[i*4+0]=3;vf->Data_index[i*4+1]=i*3+0;vf->Data_index[i*4+2]=i*3+1;vf->Data_index[i*4+3]=i*3+2;
        i++;
    }
    vf->color_rows=vf->Data_index_rows;
    vf->set_color(vf,0.9,0.0,0.0,1.0);

    vf->normal_rows=vf->Data_index_rows;
    
    free_node(n);
}

void test_output_points(char* filename ,char* outfile_name)
{
#ifdef _WIN32
    FILE *outfile;
    fopen_s(&outfile,outfile_name,"w");
#else
    FILE *outfile=fopen(outfile_name,"w");
#endif
    if(!outfile)
    {
        printf("can't open this file\n");
        return ;
    }
    fseek(outfile,0,SEEK_SET);
    fprintf(outfile,"OFF\n");

    Mesh mesh;
    Mesh_init(&mesh);
    _ReadOff_(&mesh,"face.off",3);

    fprintf(outfile,"%d %d %d\n",mesh.num_v(&mesh),0,0);
    for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
    {
         fprintf(outfile,"%lf %lf %lf\n",quote(vit)->point[0],quote(vit)->point[1],quote(vit)->point[2]);

    } 




    Mesh_free(&mesh);
    fclose(outfile);

}


static inline int is_file_obj(char* file)
{
    int len=strlen(file);
    if(file[len-1]=='j'&&file[len-2]=='b'&&file[len-3]=='o'&&file[len-4]=='.')
    {
        return 1;
    }  
    return 0;
}



void test1()
{
    int count=0;
    //char file_dir[]="E:/A2-1008-compressed";
    char file_dir[]="A1-1008-compressed";

    char** files=search_dir_files(file_dir,&count); 
    Node*node_triangles=NULL;
    for(int i=0;i<count;i++)
    {
        if(is_file_obj(files[i]))
        {
            Mesh mesh;
            Mesh_init(&mesh);
            printf("%s\n",files[i]); 
            char file_path[200]={0};
            strcat(file_path,file_dir);strcat(file_path,"/");
            strcat(file_path,files[i]);
            _ReadObj_(&mesh,file_path);
            for(auto cit=mesh.c_begin(&mesh);cit!=mesh.c_end(&mesh);cit++)
            {
                double* triangle=(double*)malloc(sizeof(double)*9);
                for(int i=0;i<3;i++)
                {
                    triangle[i*3+0] =quote(cit)->vertices[i]->point[0]; 
                    triangle[i*3+1] =quote(cit)->vertices[i]->point[1]; 
                    triangle[i*3+2] =quote(cit)->vertices[i]->point[2]; 
                }
                node_triangles=node_overlying(node_triangles,triangle);
            }
            Mesh_free(&mesh);
        }
    }  

    double* triangles=(double*)malloc(sizeof(double)*9*node_size(node_triangles));
    int sum=0;
    for(Node*nit= node_triangles;nit!=NULL;nit=(Node*)(nit->Next))
    {
        memmove(&triangles[sum*9],nit->value,sizeof(double)*9); 
        sum++;
    }  
    printf("len :%d\n",sum*9);
    //double ps[12]={6.42724,26.3074,-53.6318,11.044,279.27,-54.1057,150.311,261.695,-50.4808,178.763,30.4529,-53.4603}; 
    double ps[12]={-5.384397,130.199790,-50.412270,-5.352516,-91.639106,-50.094821,378.768309,-97.308387,-49.800767,378.728873,130.144239,-50.050807}; 


   // double ps[12]={6.42724,26.3074,0.0,11.044,279.27,0.0,150.311,261.695,0.0,178.763,30.4529,0.0}; 
    
    double*re= compute_volume_difference_of_two_mesh(triangles,node_size(node_triangles)*9,ps,12);





    Viewer_World_Manager vwm;
    viewer_world_manager_init(&vwm);
    Viewer_World *vw=vwm.create_world(&vwm,NULL);

    add_default_somethings(vw);  
    //test_show_mesh_triangles(vw,node_triangles);
    //test_show_triangles_set(vw,&re[2],re[1]);
    test_show_points_set(vw,&re[2],re[1]);
    
    Viewer_Opengl_Interpreter voi;
    viewer_opengl_interpreter_initn(&voi,&vwm);
    voi.interpreter(&voi);


    viewer_opengl_interpreter_free(&voi);
    viewer_world_manager_free(&vwm); 
    free(re);
    free(triangles);
    free_node_value(node_triangles);
    free_node(node_triangles);

    for(int i=0;i<count;i++)
    {
        printf("%s\n",files[i]);
        SAFE_FREE(files[i]); 
    }
    free(files);
}

int main()
{
    //test_output_points("face.off","face_points.off"); 
    //test_two_method("face.off","face_points.off");

    test1();
   
    printf("end\n"); 
    return 0;

}